//                       
//                  |           |           |           |
//                --^-----------^-----------^-----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     | (i,j+1,k) |     :     |
//                  o     >     o    w_N    o     >     o   
//                  |     :     |     :     |     :     |
//                  |     :     |     :     |     :     |
//                --^---------- 3 -- v_n -- 4 ----------^--  4 = v(i+1, j  ,k)
//                  |     :     |     :     |     :     |    3 = v(i  , j  ,k)
//                  |     :     |     :     |     :     |    2 = v(i+1, j-1,k)
//                  o    w_B   w_b   w_P   w_f   w_F    o    1 = v(i  , j-1,k)
//                  | (i,j,k-1) |  (i,j,k)  | (i,j,k+1) |
//                  |     :     |     :     |     :     |
//                --^---------- 1 -- v_s -- 2 ----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |     :     |     :     |
//                  o     >     o    w_S    o     >     o   
//                  |     :     | (i,j-1,k) |     :     |
//                  |     :     |     :     |     :     |
//                --^-----------^-----------^-----------^--
//                  |           |           |           | 
//                   
//   w_b = ( w(i,j,k-1) + w(i,j,k) ) / 2     
//   w_f = ( w(i,j,k+1) + w(i,j,k) ) / 2
//   v_n = ( v(i,j,k) + v(i+1,j,k) ) / 2   
//              3           4
//   v_s = ( v(i,j-1,k) + v(i+1,j-1,k) ) / 2
//              1           2

//   w_e = ( u(i,j,k) + u(i,j+1,k) ) / 2
//   u_w = ( u(i-1,j,k) + u(i-1,j+1,k) ) / 2
// 

namespace Tuna {

template<class Tprec, int Dim>
inline
bool CDS_ZLES<Tprec, Dim>::calcCoefficients(const ScalarField &nut) {
    Tprec dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    Tprec dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    Tprec ce, cw, cn, cs, cf, cb;
    Tprec nutinter;
    Tprec dxyz_dt = dx * dy * dz / dt;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = ( u(i,j,k) + u(i,j+1,k) ) * 0.5 * dyz;
		cw = ( u(i-1,j,k) + u(i-1,j+1,k) ) * 0.5 * dyz;
		cn = ( v(i,j,k) + v(i+1,j,k) ) * 0.5 * dxz;
		cs = ( v(i,j-1,k) + v(i+1,j-1,k) ) * 0.5 * dxz;
		cf = ( w(i,j,k) + w(i,j,k+1) ) * 0.5 * dxy;
		cb = ( w(i,j,k) + w(i,j,k-1) ) * 0.5 * dxy;

//
// nut is calculated on center of volumes, therefore, nut
// must be staggered in z direction:		
		nutinter = 0.5 * ( nut(i,j,k) + nut(i,j,k+1) );

		aE (i,j,k) = (Gamma + nutinter) * dyz_dx - ce * 0.5;
		aW (i,j,k) = (Gamma + nutinter) * dyz_dx + cw * 0.5;
		aN (i,j,k) = (Gamma + nutinter) * dxz_dy - cn * 0.5;
		aS (i,j,k) = (Gamma + nutinter) * dxz_dy + cs * 0.5;
		aF (i,j,k) = 2 * (Gamma + nutinter) * dxy_dz - cf * 0.5;
		aB (i,j,k) = 2 * (Gamma + nutinter) * dxy_dz + cb * 0.5;

		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) +
		             aN (i,j,k) + aS (i,j,k) +
		             dxyz_dt;
//		aP (i,j,k) /= alpha;  // under-relaxation
//		+ (ce - cw) + (cn -cs) + (cf - cb);	    
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.
		sp(i,j,k) = w(i,j,k) * dxyz_dt - 
		    ( p(i,j,k+1)- p(i,j,k) ) * dxy + 
		    nutinter * ( (u(i,j,k+1) - u(i,j,k) - 
				  u(i-1,j,k+1) + u(i-1,j,k)) * dy +
				 (v(i,j,k+1) - v(i,j,k) - 
				  v(i,j-1,k+1) + v(i,j-1,k)) * dx );
		
//		    w(i,j,k) * (1-alpha) * aP(i,j,k)/alpha;// under-relaxation
	}    

    calc_dw_3D();
    applyBoundaryConditions3D();
  
    return 1;
}

} // Tuna namespace














